#setwd('/data1/bsi/bioinf_int/s106381.borawork/data_HTHGU/HTHGU_GeneSumm/samplefiles_out')
args <- commandArgs(TRUE)
setwd(args[1])
library(genefilter, lib.loc=args[2])
library(gcrma, lib.loc=args[2])
library(Biobase, lib.loc=args[2])
library(preprocessCore, lib.loc=args[2])
library(MASS, lib.loc=args[2])
library(affy, lib.loc=args[2])
library(simpleaffy, lib.loc=args[2])
library(frma, lib.loc=args[2])
library(hgu133plus2geneccdscdf, lib.loc=args[2])
library(hgu133plus2frmavecs,lib.loc=args[2])
#library(hgu133afrmavecs, lib.loc=args[2])
#library(equalprobesorihgu133ageneccdscdf, lib.loc=args[2])
library(affyPLM, lib.loc=args[2])
data <- ReadAffy()
data@cdfName<-"hgu133plus2geneccds"
#data@cdfName<-"equalprobesorihgu133ageneccds"
#script start
#######################################################################################

frma1 <- function(object, background="rma", normalize="quantile", summarize="robust_weighted_average", input.vecs=list(normVec=NULL, probeVec=NULL, probeVarBetween=NULL, probeVarWithin=NULL, probesetSD=NULL), output.param=NULL, verbose=FALSE){

  if(!is(object, "AffyBatch")) stop(paste("argument is", class(object), "frma requires AffyBatch"))

  if(!background %in% c("none", "rma")) stop("background must be either none or rma")
  if(!normalize %in% c("none", "quantile")) stop("normalize must be either none or quantile")
  if(!summarize %in% c("median_polish", "average", "median", "weighted_average", "robust_weighted_average", "batch")) stop("summarize must be one of: median_polish, average, median, weighted_average, robust_weighted_average, batch")

  if(summarize=="batch" & (any(input.vecs$probeVarBetween==0) | any(input.vecs$probeVarWithin==0))) stop("If summarize method is batch then probeVarBetween and probeVarWithin must be greater than zero for all probes.")
  
  cdfname <- cleancdfname(cdfName(object))
  platform <- gsub("cdf","",cdfname)

  if(summarize == "median_polish") output <- frmaMedPol(object, background, normalize, input.vecs, verbose)
  if(summarize %in% c("average", "median", "weighted_average", "robust_weighted_average")) output <- frmaRobReg1(object, background, normalize, summarize, input.vecs, output.param, verbose)
  if(summarize == "batch") output <- frmaBatch1(object, background, normalize, input.vecs, output.param, verbose)

  if("stderr" %in% output.param) stderr <- output$stderr else stderr <- matrix(nrow=0, ncol=0)
  if("weights" %in% output.param) w <- output$weights else w <- matrix(nrow=0, ncol=0)
  if("residuals" %in% output.param) r <- output$residuals else r <- matrix(nrow=0, ncol=0)

  R.model <- PLM.designmatrix3(object)
  
  new("PLMset",
      chip.coefs=output$exprs,
      weights=list("PM.weights"=w, "MM.weights"=matrix(nrow=0, ncol=0)),
      se.chip.coefs=stderr,
      residuals=list("PM.resid"=r, "MM.resid"=matrix(nrow=0, ncol=0)),
      cdfName = object@cdfName,
      phenoData = phenoData(object),
      annotation = annotation(object),
      experimentData = experimentData(object),
      nrow= object@nrow,
      ncol= object@ncol,
      narrays=length(object),
      model.description = list("R.model"=R.model))
}
########################################################################################
frmaRobReg1 <- function(object, background, normalize, summarize, input.vecs, output.param, verbose){

  cdfname <- cleancdfname(cdfName(object))
  cdfname<-"hgu133plus2"
 
platform <- gsub("cdf","",cdfname)
  
  if(background == "rma"){
    if(verbose) message("Background Correcting ...\n")
    object <- bg.correct.rma(object)
  }
  
  if(is.null(input.vecs$normVec) | is.null(input.vecs$probeVec) | is.null(input.vecs$probeVarWithin) | is.null(input.vecs$probeVarBetween) | (summarize=="robust_weighted_average" & is.null(input.vecs$probesetSD))){
    pkg <- paste(platform, "frmavecs", sep="")
    require(pkg, character.only=TRUE, quiet=TRUE) || stop(paste(pkg, "package must be installed first"))
    data(list=eval(pkg))

    if(is.null(input.vecs$normVec)) input.vecs$normVec <- get(pkg)$normVec
    if(is.null(input.vecs$probeVec)) input.vecs$probeVec <- get(pkg)$probeVec
    if(is.null(input.vecs$probeVarWithin)) input.vecs$probeVarWithin <- get(pkg)$probeVarWithin
    if(is.null(input.vecs$probeVarBetween)) input.vecs$probeVarBetween <- get(pkg)$probeVarBetween
    if(is.null(input.vecs$probesetSD) & summarize=="robust_weighted_average") input.vecs$probesetSD <- get(pkg)$probesetSD
  }

  if(normalize == "quantile"){
    if(verbose) message("Normalizing ...\n")
    pm(object) <- normalize.quantiles.use.target(pm(object), input.vecs$normVec)
  }

  if(verbose) message("Summarizing ...\n")
  pns <- probeNames(object)
  pms <- log2(pm(object))
  
  if(summarize == "average"){  
    exprs <- subColSummarizeAvg((pms-input.vecs$probeVec), pns)
    weights <- NULL
    stderr <- NULL
  }

  if(summarize == "median"){  
    exprs <- subColSummarizeMedian((pms-input.vecs$probeVec), pns)
    weights <- NULL
    stderr <- NULL
  }
  
  if(summarize == "weighted_average"){  
    w <- 1/(input.vecs$probeVarWithin + input.vecs$probeVarBetween)
    if(any(input.vecs$probeVarWithin==0) | any(input.vecs$probeVarBetween==0)) warning("Either probeVarWithin or probeVarBetween is 0 for some probes -- setting corresponding weights to 1")
    w[w==Inf] <- 1
    exprs <- subColSummarizeAvg((pms-input.vecs$probeVec)*w, pns)
    W <- as.vector(rowsum(w, pns, reorder=FALSE))
    exprs <- (exprs/W)*as.vector(rowsum(rep(1,length(pns)),pns,reorder=FALSE))
    weights <- NULL
    stderr <- NULL
  }

  if(summarize == "robust_weighted_average"){
    w <- 1/(input.vecs$probeVarWithin + input.vecs$probeVarBetween)
    if(any(input.vecs$probeVarWithin==0) | any(input.vecs$probeVarBetween==0)) warning("Either probeVarWithin or probeVarBetween is 0 for some probes -- setting corresponding weights to 1")
    w[w==Inf] <- 1
    N <- 1:dim(pms)[1]
    S <- split(N, pns)
    fit <- lapply(1:length(S), function(i) {
	s <- S[[i]]
	rwaFit2(pms[s,, drop=FALSE], w[s], input.vecs$probeVec[s], input.vecs$probesetSD[i])
    })
    names(fit) <- unique(pns)
    exprs <- matrix(unlist(lapply(fit, function(x) x$Estimates)), ncol=ncol(pms), byrow=TRUE)
    rownames(exprs) <- names(fit)
    colnames(exprs) <- colnames(pms)
    if("weights" %in% output.param){
      weights <- matrix(unlist(lapply(fit, function(x) x$Weights)), ncol=ncol(pms), byrow=TRUE)
      rownames(weights) <- pns
      colnames(weights) <- sampleNames(object)
    } else weights <- NULL
    if("stderr" %in% output.param){
      stderr <- matrix(unlist(lapply(fit, function(x) x$StdErrors)), ncol=ncol(pms), byrow=TRUE)
      rownames(stderr) <- names(fit)
      colnames(stderr) <- sampleNames(object)
    } else stderr <- NULL
  }

  if("residuals" %in% output.param){
    residuals <- apply(exprs,2,function(x) rep(x, table(pns)))
    residuals <- (pms-input.vecs$probeVec) - residuals
    rownames(residuals) <- pns
    colnames(residuals) <- sampleNames(object)
  } else residuals <- NULL
  
  colnames(exprs) <- sampleNames(object)
  
  return(list(exprs=exprs, stderr=stderr, weights=weights, residuals=residuals))
}

rwaFit2 <- function(x1, x2, x3, x4){
  ncols <- ncol(x1)
  w.tmp <- x2/max(x2)
  w.tmp <- matrix(rep(w.tmp, ncols), ncol=ncols)
  pe.tmp <- x3
  pe.tmp[1] <- pe.tmp[1]-sum(pe.tmp)
  rcModelWPLM(y=x1, w=w.tmp, row.effects=pe.tmp, input.scale=x4)
}


#########################################################################################
frmaBatch1 <- function(object, background, normalize, input.vecs, output.param, verbose){

  cdfname <- cleancdfname(cdfName(object))
  cdfname<-"hgu133plus2"
	platform <- gsub("cdf","",cdfname)
  
  if(background == "rma"){
    if(verbose) message("Background Correcting ...\n")
    object <- bg.correct.rma(object)
  }

  if(is.null(input.vecs$normVec) | is.null(input.vecs$probeVec) | is.null(input.vecs$probeVarWithin) | is.null(input.vecs$probeVarBetween)){
    pkg <- paste(platform, "frmavecs", sep="")
    require(pkg, character.only=TRUE, quietly=TRUE) || stop(paste(pkg, "package must be installed first"))
    data(list=eval(pkg))

    if(is.null(input.vecs$normVec)) input.vecs$normVec <- get(pkg)$normVec
    if(is.null(input.vecs$probeVec)) input.vecs$probeVec <- get(pkg)$probeVec
    if(is.null(input.vecs$probeVarWithin)) input.vecs$probeVarWithin <- get(pkg)$probeVarWithin
    if(is.null(input.vecs$probeVarBetween)) input.vecs$probeVarBetween <- get(pkg)$probeVarBetween
  }

  if(normalize == "quantile"){
    if(verbose) message("Normalizing ...\n")
    pm(object) <- normalize.quantiles.use.target(pm(object), input.vecs$normVec)
  }

  if(verbose) message("Summarizing ...\n")
  pns <- probeNames(object)
  pms <- log2(pm(object))
  
  tmp <- split(data.frame(pms, input.vecs$probeVec, input.vecs$probeVarWithin, input.vecs$probeVarBetween), pns)
  fit <- lapply(tmp, batchFit)
  
  exprs <- matrix(unlist(lapply(fit, function(x) x$exprs)), ncol=ncol(pms), byrow=TRUE)
  rownames(exprs) <- names(fit)
  colnames(exprs) <- colnames(pms)

  if("weights" %in% output.param){
    weights <- matrix(unlist(lapply(fit, function(x) x$w)), ncol=ncol(pms), byrow=TRUE)
    rownames(weights) <- pns
    colnames(weights) <- sampleNames(object)
  } else weights <- NULL

  if("stderr" %in% output.param){
    stderr <- unlist(lapply(fit, function(x) x$se))
    names(stderr) <- names(fit)
  } else stderr <- NULL

  if("residuals" %in% output.param){
    residuals <- apply(exprs,2,function(x) rep(x, table(pns)))
    residuals <- (pms-input.vecs$probeVec) - residuals
    rownames(residuals) <- pns
    colnames(residuals) <- sampleNames(object)
  } else residuals <- NULL
  
  colnames(exprs) <- sampleNames(object)
  
  return(list(exprs=exprs, stderr=stderr, weights=weights, residuals=residuals))
}




rcModelWPLM<-function (y, w, row.effects = NULL, input.scale = NULL)
{
    if (!is.matrix(y))
        stop("argument should be matrix")
    if (is.vector(w)) {
        if (length(w) != prod(dim(y))) {
            stop("weights are not correct length")
        }
    }
    else if (is.matrix(w)) {
        if (!all(dim(w) == dim(y))) {
            stop("weights should be same dimension as input matrix")
        }
    }
    if (any(w < 0)) {
        stop("weights should be no negative")
    }
    PsiCode <- 0
    PsiK <- 1.345
    if (is.null(row.effects)) {
        .Call("R_wrlm_rma_default_model", y, PsiCode, PsiK, as.double(w),
            input.scale, PACKAGE = "preprocessCore")
    }
    else {
        if (length(row.effects) != nrow(y)) {
            stop("row.effects parameter should be same length as number of rows")
        }
        if (abs(sum(row.effects)) > 143 * .Machine$double.eps) {
            stop("row.effects should sum to zero")
        }
        .Call("R_wrlm_rma_given_probe_effects", y, as.double(row.effects),
            PsiCode, PsiK, as.double(w), input.scale, PACKAGE = "preprocessCore")
    }
}



#script end

affy.dc5 <- mas5calls(data)
affy.pval <- assayData(affy.dc5)[["se.exprs"]]
#pval <- detection.p.val(data)
#pval.c <- pval$call
#pval.p <- pval$pval
#pval.pc <- cbind(pval.p,pval.c)
object <- frma1(data, summarize = "robust_weighted_average")
mydata <- as.vector(strsplit(args[1], '/')[[1]])
output.file=paste(mydata[length(mydata)],".normalized",sep="")
output.file.call = paste(output.file,".call",sep="")
#write.exprs(object,output.file,sep="\t")
write.table(coefs(object),output.file,sep="\t",quote=FALSE)
write.table(affy.pval,output.file.call,sep="\t",quote=FALSE)
#dat<-read.table(output.file,sep="\t",header=T)
#final_dat<-cbind(dat,pval.pc)
#nam <-as.vector(colnames(final_dat))
#nam[1]<-"Probeset"
#nam[5]<-"Sample Name"
#nam[6]<-"Affy_Platform"
#final_dat<-cbind(final_dat,nam[2],annotation(data))
#colnames(final_dat)<-nam
#write.table(final_dat,output.file,sep="\t",row.names=F,quote=F)
q()
n
